*                           R E S A M P L E . F T N
*                           ***********************

*     NOTE: REMOVE COMMENTS BEFORE RUNNING.  DATA ARRAYS NEED THE SPACE.

*     For purposes of time series analysis, the number  of  data  points
*     must be an exact power of 2, at least when Fast Fourier Transform,
*     and related filter algorithms are used.  However,  a  time  series
*     will generally not consist of exactly the correct number  of  data
*     points, particularly if the data has been obtained  through  tests
*     or experimentation.  It is common to simply truncate the data file
*     or pad zeros to the end in order to obtain the required number  of
*     points.  Neither of these methods is entirely satisfactory.

*     The present program accepts a time series data file as  input  and
*     creates from it  a  new  file  containing  any  number  of  points
*     desired, including powers of  2.   When  an  increased  number  of
*     points is requested, direct spline interpolation is used to obtain
*     the new data points which will lie between the  old  data  points.
*     When a smaller number of data points is requested, The time series
*     is first passed through a low pass filter before  resampling  with
*     the spline interpolation.  The filter cut-off frequency is set  at
*     less than half the new sample  frequency  to  avoid  any  aliasing
*     difficulties.

*     The program begins by requesting the name of the source data file.
*     It is assumed that the source file will have been formatted by the
*     program 'dataform.ftn' or 'datout'.  The four  header  records  of
*     this file will be read and displayed to show the  number  of  data
*     points, the time origin and the sample interval.  A  request  will
*     then be made for  the  number  of  data  points  desired  and  the
*     destination  file  name.   The  new  sample  interval  and  cutoff
*     frequency will be displayed and the processing will proceed.   The
*     new sample interval dt' is calculated from the old sample interval'
*     dt using the relation:

*          dt' = dt * ( n - 1 ) / ( n' - 1 )

*     where n is the number of data points in the source file and n'  is'
*     the number of points in the new file.  The  cutoff  frequency   fc
*     for the original data file is calculated by assuming it  is  equal
*     to half the sample frequency.  Thus:

*          fc =  (1/2) * 1 / dt

*     and this is expressed in cycles per second if dt  is  measured  in
*     seconds.  if n' is less than n,  a  low  pass  filter  is  applied'
*     before interpolation.  The new cutoff frequency will then be:

*          fc' = (1/2) * 1 / dt'

*     When this is expressed in  the  usual  normalized  form  with  the
*     sample interval taken as unity, and  the  frequency  expressed  in
*     radians per sample interval then:

*          omega' = fc' * pi * dt

*                 = (1/2) * pi * ( n' - 1 ) / ( n - 1)'

*     P R O G R A M   S T R U C T U R E
*     *********************************

*     The program roughly conforms to the following steps:

*          Request name of source data file
*          open Source data file
*          read source data file header
*          Display header information
*          read source data file into data array
*          close source file
*          Request destination file information
*          Check that n' .gt. 2  and n' .ne. n

*          if  n' .LT. n  THEN'
*               Calculate new file cutoff frequency and filter length
*               Display destination file header information
*               Extend data array at both ends
*               Apply lowpass filter to data array
*          else
*               Display destination file information
*          endif

*          Save the filtered data in a temporary file
*          execute subroutine to calculate spline sigma array
*          Re-read the filtered data from the temporary file

*          open output file
*          write output file header records
*          Initialize data record counter
*          while  record counter .lt. next to last record
*               Calculate interpolated data one record at a time
*               write output record
*          endwhile
*          Pad last record with zeros if necessary
*          write last Record
*          close output file

*          end


*                        L O W P A S S    F I L T E R
*                        ****************************

*     This subroutine replaces the  series   x   of  length   n    by  a
*     smoothed version obtained with a  low-pass  digital  filter.   The
*     parameter  m  is the maximum length of  x , that is, the dimension
*     assigned to  x  in the calling  program.   The  filter  length  is
*     (2 * ns + 1)  terms and the filter  weights  represent  the  least
*     squares approximation to the cutoff filter with  cutoff  frequency
*     w .  The cutoff frequency  w  is expressed in radians per  sample.
*     The filter coefficients are proportional to:

*               sin( i * w )    sin( 2 * pi * i / ( 2 * ns + 1) )
*               -----------  *  --------------------------------
*                   i * w          2 * pi * i / ( 2 * ns + 1 )

*     where the first term  is  simply  the  Fourier  Transform  of  the
*     rectangular pass band with cutoff frequency  w ,  and  the  second
*     term is a convergence factor to  limit  the  Gibb's  oscillations.'
*     The coefficients are  normalized  so  that  their  sum  is  unity,
*     providing unity gain for a constant input array, that is, at  zero
*     frequency.  The transfer function for  this  filter  will  have  a
*     transition band of frequency width given by  4 * pi / (2 * ns + 1)
*     which suggests that a  large  value  of   ns   would  reduce  this
*     transition band width.

*     to avoid loss of data due to the finite length of the filter,  the
*     input  array  is  extended  symmetrically  at  both  ends  by   ns
*     positions.  This approach is not perfect since it  will  introduce
*     some phase distortion in the filtered results at both  ends.   The
*     extended array occupies  n + 2 * ns   positions  and  the  maximum
*     size allowed for  x  must be at least as  large  as  the  extended
*     array.  The filtered result is left in the  lower   n   positions.
*     The size of  ns  is a balance  between  the  desire  for  a  small
*     frequency transition band and a desire  to  limit  the  distortion
*     length  at  the  array  ends  as  well  as   to   accelerate   the
*     computations.  A cost may be defined consisting of the sum of  the
*     proportion of the frequency transition band to the pass band,  and
*     the proportion of the distortion length to the array  length.   An
*     "optimum" value of  ns  may then be obtained  by  minimizing  this
*     cost.  Thus:

*                         4 * pi / ( 2 * ns + 1 )           2 * ns
*         cost   =   A *  -----------------------   +   B * ------
*                                  w                          n

*     where  A  and  B  are weights which may  be  used  to  adjust  the
*     importance of each term.  The minimization results in  the  filter
*     half length of:

*          ns    =    ( SQRT( (A / B) * 4 * pi * n / w ) - 1 ) / 2

*     and the smallest integer part of this expression is taken.

*     The particular filter chosen  here  is  a  two  sided,  symmetric,
*     non-recursive filter.  This  ensures  that  phase  distortion  and
*     delays will not be introduced, except  at  the  end  points.   The
*     filtering could be accomplished more rapidly  if  a  single  sided
*     non-recursive filter were employed since FFT methods could then be
*     used.  However, the time delays introduced this way could make  it
*     difficult to relate the filtered data to real time.   Consequently
*     the two sided filter has been chosen.


*                  S P L I N E    I N T E R P O L A T I O N
*                  ****************************************

*     The spline  interpolation  subroutine  is  based  on  the  program
*     published  in  the  book  "Computer   Methods   for   Mathematical"
*     Computations" by Forsythe, Malcolm and Moler, Prentice-Hall  1977,"
*     Page  76.   The  program  has  been  modified  to  reduce  storage
*     requirements  by  capitalizing  on  the  fact  that  the  data  is
*     available at equally spaced intervals.  Further more, rather  than
*     calculating and storing the three cubic coefficients for each data
*     point, only the so called 'sigma' array is stored.  That  is,  the
*     interpolated points are obtained from the following  interpolation
*     equation:

*         x(z)  =  w * x(i+1)  +  w' * x(i)'

*                    +  h*h * [w*(w*w-1)*s(i+1) + w'*(w'*w'-1)*s(i)]'

*     where
*                 h  = sample interval

*                 w  = z - h * i

*                 w' = 1 - w'

*               s(i) = i'th sigma quantity'

* ************************************************************************

      real x(1700), s(1700), tzero, dt, t(5)
      character source$file, destination$file, file$id, var$fmt, garbag

* DIMENSIONS OF  x  AND  s  AND THE PARAMETER  m   DETERMINE MAX ARRAY SIZE.

      m = 1700
      pi = 4.0 * atan(1.0)

* CLEAR SCREEN, IDENTIFY PROGRAM AND REQUEST DATA FILE NAME.

      write(6,"'1'")
      print,"        P R O G R A M    R E S A M P L E . F T N"
      print,"        ****************************************"
      print
      print," Enter the source file name"
      read, source$file

* READ THE INPUT DATA FILE.

      call data$input(source$file, file$id, m, var$fmt, tzero, dt, fc, x, n)

* HOW MANY POINTS IN RE-SAMPLED FILE?  ENSURE NUMBER IS IN ACCEPTABLE RANGE.

      loop
          write(6,1) char(11),char(11)
    1     format(1x,2a1,"Enter number of data points for the destination file")
          read, nprime
          quitif (nprime .gt. 2) .and. (nprime .ne. n)
          write(6,"'+          '")
      endloop

* ENTER AND DISPLAY OUTPUT FILE HEADER INFORMATION.
* APPLY LOW PASS FILTER IF RESAMPLING AT A REDUCED RATE.

      print,"Enter destination file name"
      read, destination$file
      file$id = file$id//"- RESAMPLED"
      dtprime = dt * float(n - 1) / float(nprime - 1)
      print
      print,"DESTINATION FILE PARAMETERS"
      if nprime .lt. n
          fcprime = 3600.0 * 0.5 / dtprime
          call display(file$id, nprime, var$fmt, tzero, dtprime, fcprime)
          w = 0.5 * pi * float(nprime-1) / float(n-1)
          call lowpass(x, n, m, w)
      else
          fcprime = fc
          call display(file$id, nprime, var$fmt, tzero, dtprime, fcprime)
      endif

* CALCULATE SPLINE INTERPOLATION ARRAY FOR EQUALLY SPACED DATA BUT FIRST
* SAVE THE FILTERED  x  ARRAY IN A TEMPORARY FILE SINCE splineq WILL
* DESTROY IT.  RE-READ THE DATA WHEN splineq IS FINISHED.

      open(unit=8,file="temporary.dat")
      write(8,var$fmt) (x(i),i=1,n)
      close(unit=8)
      call splineq(n, m, x, s)
      open(unit=8,file="temporary.dat")
      read(8,var$fmt) (x(i),i=1,n)
      close(unit=8)

* WRITE THE FOUR OUTPUT FILE HEADER RECORDS.  USE THE STANDARD FORMATS.

      open(unit=8,file=destination$file)
      write(8,"a79,'1'") file$id
      write(8,"i5,74x,'2'") nprime
      write(8,"a79,'3'") var$fmt
      write(8,"3f10.5,49x,'4'") tzero, dtprime, fcprime

* OBTAIN INTERPOLATED POINTS AND WRITE TO FILE, RECORD BY RECORD.

      l = 4
      factor = float(n-1)/float(nprime-1)
      jend = nprime/5
      do j = 1 , jend
          k = 5 * (j - 1) - 1
          do i = 1 , 5
              u = 1.0 + float(k+i) * factor
              t(i) = val(n, u, x, s)
          enddo
          l = j + 4
          write(8,"5e15.8,i5") (t(i),i=1,5), l
      enddo

* PAD LAST RECORD WITH ZEROS IF NEEDED.

      k = 5 * jend
      if k .lt. nprime
          k = k - 1
          do i = 1 , 5
              t(i) = 0.0
              u = 1.0 + float(k+i) * factor
              if (k+i) .lt. nprime
                  t(i) = val(n, u, x, s)
              endif
          enddo
          l = l + 1
          write(8,"5e15.8,i5") (t(i),i=1,5), l
      endif

      close(unit=8)
      end

      subroutine data$input(file$name, file$id, m, var$fmt, tzero, dt, fc, x, n)
      character file$name, file$id, var$fmt, two$characters, two$blanks
      integer n, character$count
      real tzero, dt, fc, x(m)
      open(unit=8,file=file$name)

* READ FILE ID AND REMOVE TRAILING BLANKS.

      read(8,"a72") file$id
      character$count = 0
      two$blanks = "  "
      while character$count .lt. 72
          character$count = character$count + 1
          two$characters = substr(file$id, character$count, 2)
          quitif two$characters = two$blanks
      endwhile
      file$id = substr(file$id, 1, character$count)

* READ REMAINING THREE HEADER RECORDS.

      read(8,"i5") n
      read(8,"a9") var$fmt
      read(8,"3f10.5") tzero, dt, fc

* CALCULATE ASSUMED FILE CUTOFF FREQUENCY IF NOT RECORDED.

      if  fc = 0.0
          fc = 3600.0 * 0.5 / dt
      endif

* DISPLAY THE HEADER INFORMATION.

      call display(file$id, n, var$fmt, tzero, dt, fc)

* READ THE DATA.

      read(8,var$fmt) (x(j),j=1,n)
      close(unit=8)
      return
      end


      subroutine display(id,n,vform,t,dt,f)
      character id,vform
      print
      print,"      File Identifier: ", id
      print,"Number of Data Points: ", n
      print,"          Data Format: ", vform
      print,"          Time Origin: ", t
      print,"            Time Step: ", dt
      write(6,1) f
    1 format("      Cutoff frequency: ", f6.2, " Cycles per Hour"///)
      return
      end


      subroutine lowpass(x, n, m, w)
      real x(m), h(100)
      integer hsize
      hsize = 100
      pi = 4.0 * atan(1.0)

* CALCULATE FILTER HALF LENGTH.

      A = 1.0
      B = 10.0
      ns = ( sqrt( (A / B) * 4 * pi * float(n) / w ) -1 ) / 2.0
      length = 2 * ns + 1
      if  length .gt. hsize
          print,"Filter length =",length," exceeds dimension specification"
          stop
      endif
      if  (n + 2 * ns) .gt. m
          print,"Extended  x  array exceeds available dimension"
          stop
      endif

*  EXTEND THE ARRAY.

      n1 = n+ns+1
      n2 = n+1
      do i = 1 , n
          x(n1-i) = x(n2-i)
      enddo
      n1 = ns+1
      n2 = n+ns
      do i = 1 , ns
          x(n1-i) = x(n1+i)
          x(n2+i) = x(n2-i)
      enddo

* CALCULATE FILTER WEIGHTS.

      hzero = w / pi
      con = 2.0 * pi / float(2*ns+1)
      sum = hzero
      do i = 1 , ns
          h(i) = hzero * sinc( float(i) * w ) * sinc( float(i) * con )
          sum  = sum + 2.0 * h(i)
      enddo
      hzero = hzero / sum
      do i = 1 , ns
          h(i) = h(i) / sum
      enddo

* APPLY FILTER WEIGHTS.

      do i = 1 , n
          ins  = i + ns
          temp = hzero * x(ins)
          do j = 1 , ns
              temp = temp + (x(ins+j) + x(ins-j)) * h(j)
          enddo
          x(i) = temp
      enddo
      return
      end

      function sinc(z)
          sinc = sin(z) / z
      end



      subroutine splineq(n, m, x, s)
      real x(m), s(m)

* NOTE: x  ARRAY IS DESTROYED BY THIS ROUTINE.

      if (n .le. 3)
          print," Array length is less than or equal to 3.  This"
          print," spline routine needs at least 4 points."
          print," Execution has been terminated."
          stop
      endif

*  SETUP TRIDIAGONAL SYSTEM FOR EQUALLY SPACED INTERVALS.

      nm1 = n - 1
      s(2) = x(2)-x(1)
      bi = 4.0
      do i = 2 , nm1
          s(i+1) = x(i+1) - x(i)
          s(i) = s(i+1) - s(i)
      enddo

*  END CONDITIONS.

      b1 = -1.0
      bn = -1.0
      s(1) = ( s(3)   - s(2)   ) / 2.0
      s(n) = ( s(n-1) - s(n-2) ) / 2.0
      s(1) =  s(1) / 3.0
      s(n) = -s(n) / 3.0

*  FORWARD ELIMINATION.  THIS IS WHERE THE  x  ARRAY IS USED FOR STORAGE
*  OF THE DIAGONAL  b  ELEMENTS.

      b = b1
      x(1) = b
      do i = 2 , n
          t = 1.0/b
          b = bi - t
          x(i) = b
          s(i) = s(i)-t*s(i-1)
      enddo
      bn = bn - t

*  BACK SUBSTITUTION.

      s(n) = s(n)/bn
      b = 1.0/t
      do j = 1 , nm1
          i = n-j
          b = x(i)
          s(i) = (s(i)-s(i+1))/b
      enddo
      return
      end

      real function val(n, u, x, s)
      real x(n), s(n), u
      integer i
      i = u
      if (i .ge. n) i = n - 1
      w = u - float(i)
      wbar = 1.0 - w
      temp = w * x(i+1) + wbar * x(i)
      temp = temp + w*(w+1.0)*(w-1.0)*s(i+1)
      temp = temp + wbar*(wbar+1.0)*(wbar-1.0)*s(i)
      val  = temp
      return
      end
